﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary.FFTPack
{
    /// <summary>
    /// Cosine FFT transform of a real even sequence. 
    /// </summary>
    [Serializable]
    public class RealDoubleFFTEven : RealDoubleFFTMixed
    {
        /// <summary>
        /// This field holds the norm factor.
        /// </summary>
        private double normFactor;

        /// <summary>
        /// This field holds the wave table.
        /// </summary>
        private readonly double[] wavetable;

        /// <summary>
        /// This field holds the dimension.
        /// </summary>
        private readonly int ndim;

        /// <summary>
        /// Initializes a new instance of the <see cref="RealDoubleFFTEven"/> class.
        /// </summary>
        /// <param name="n">The size of a real data sequence. When n is a multiplication of 
        /// small numbers (4, 2, 3, 5), this FFT transform is very efficient.</param>
        public RealDoubleFFTEven(int n)
        {
            ndim = n;
            normFactor = 2 * (n - 1);
            if (wavetable == null || wavetable.Length != (3 * ndim + 15))
            {
                wavetable = new double[3 * ndim + 15];
            }
            Costi(ndim, wavetable);
        }

        /// <summary>
        /// This field can be used to normalize this FFT transform. This is because a call of forward 
        /// transform followed by a call of backward transform will multiply the input sequence by the 
        /// norm factor.
        /// </summary>
        /// <value>The norm factor to use for the transform.</value>
        public double NormFactor
        {
            get { return normFactor; }
            set { normFactor = value; }
        }

        /// <summary>
        /// Forward cosine FFT transform. It computes the discrete sine transform of an odd sequence.
        /// </summary>
        /// <param name="x">An array which contains the sequence to be transformed. After FFT, x 
        /// contains the transform coeffients.</param>
        public void ForwardTransform(double[] x)
        {
            Cost(ndim, x, wavetable);
        }

        /// <summary>
        /// Backward cosine FFT transform. It is the unnormalized inverse transform of ForwardTransform().
        /// </summary>
        /// <param name="x">An array which contains the sequence to be transformed. After FFT, x 
        /// contains the transform coeffients.</param>
        public void BackTransform(double[] x)
        {
            Cost(ndim, x, wavetable);
        }

        private void Cost(int n, double[] x, double[] wtable)
        {
            int modn, i, k;
            double c1, t1, t2;
            int kc;
            double xi;
            int nm1;
            double x1h;
            int ns2;
            double tx2, x1p3, xim2;

            nm1 = n - 1;
            ns2 = n / 2;
            if (n - 2 < 0)
            {
                return;
            }
            if (n == 2)
            {
                x1h = x[0] + x[1];
                x[1] = x[0] - x[1];
                x[0] = x1h;
            }
            else if (n == 3)
            {
                x1p3 = x[0] + x[2];
                tx2 = x[1] + x[1];
                x[1] = x[0] - x[2];
                x[0] = x1p3 + tx2;
                x[2] = x1p3 - tx2;
            }
            else
            {
                c1 = x[0] - x[n - 1];
                x[0] += x[n - 1];
                for (k = 1; k < ns2; k++)
                {
                    kc = nm1 - k;
                    t1 = x[k] + x[kc];
                    t2 = x[k] - x[kc];
                    c1 += wtable[kc] * t2;
                    t2 = wtable[k] * t2;
                    x[k] = t1 - t2;
                    x[kc] = t1 + t2;
                }
                modn = n % 2;
                if (modn != 0)
                {
                    x[ns2] += x[ns2];
                }
                Rfftf1(nm1, x, wtable, n);
                xim2 = x[1];
                x[1] = c1;
                for (i = 3; i < n; i += 2)
                {
                    xi = x[i];
                    x[i] = x[i - 2] - x[i - 1];
                    x[i - 1] = xim2;
                    xim2 = xi;
                }
                if (modn != 0)
                {
                    x[n - 1] = xim2;
                }
            }
        }

        private void Costi(int n, double[] wtable)
        {
            double pi = Math.PI;
            int k, kc, ns2;
            double dt;

            if (n <= 3)
            {
                return;
            }
            ns2 = n / 2;
            dt = pi / (double) (n - 1);
            for (k = 1; k < ns2; k++)
            {
                kc = n - k - 1;
                wtable[k] = 2 * Math.Sin(k * dt);
                wtable[kc] = 2 * Math.Cos(k * dt);
            }
            Rffti1(n - 1, wtable, n);
        }
    }
}